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ABSTRACT 


Geographic  location  of  radar  emitters  is  the  process  of  estimating  an  emitter’s 
location  upon  the  surface  of  the  earth  from  direction  of  arrival  (DOA)  data  for  the 
targeted  emitter.  The  current  Emitter  Location  (EMLOC)  algorithm  utilized  by  the 
Grumman  EA-6B  Prowler  is  based  on  a  thesis  presented  by  Mr.  Richard  Opperman  in 
June  1982.  With  the  advent  of  increased  processing  demands  on  the  AN/AYK-14 
Tactical  Computer  as  part  of  recent  software  upgrades  to  the  AN/ALQ-99  Tactical 
Jamming  System,  it  was  hoped  that  a  Kalman  Filter,  or  Extended  Kalman  Filter  based 
algorithm  would  reduce  the  processing  time  and  memory  requirements  for  the  EMLOC 
algorithm.  This  thesis  compares  the  current  algorithm,  and  the  Kalman/Extended  Kalman 
Filters  in  a  tactical  scenario  to  determine  if  a  change  in  the  current  Onboard  Flight 
Program  (OFP)  should  be  recommended. 
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L  BACKGROUND 


A.  PURPOSE 

This  thesis  examines  current  and  proposed  methods  of  geographic  emitter  location 
(EMLOC)  with  respect  to  the  EA-6B  aircraft  via  the  AN/ALQ-99.  The  current  algorithm 
will  be  compared  to  a  proposed  Kalman  Filter  based  design.  The  various  alternatives  to 
EMLOC,  along  with  their  underlying  assumptions,  will  be  verified  mathematically  and  a 
proposed  test  developed. 

B.  AIRCRAFT 

The  EA-6B  aircraft  is  currently  flown  with  three  variants  of  the  Improved 
Capabilities  II  (ICAP  II)  software/hardware  upgrades.  These  are  Blocks  82,  86,  and  89A. 
The  difference  in  the  various  blocks  will  be  discussed,  however,  the  focus  of  this  paper 
will  be  Block  86/89A  configured  aircraft. 

C.  AN/ALQ-99 

The  AN/ALQ-99  is  composed  of  two  major  sections.  The  Tactical  Jamming 
System  (TJS)  carries  out  the  majority  of  the  Electronic  Attack  (EA)  capability  of  the 
system.  The  Onboard  System  (OBS)  is  effectively  the  Electronic  Support  (ES)  portion  of 
the  system.  By  utilizing  the  superhetrodyned  receivers,  the  operator  may  scan  various 
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prioritized  portions  of  the  electromagnetic  spectrum,  or  the  entire  effective  frequency 
range  for  the  receivers.  For  Block  82  aircraft,  as  receivers  are  scanned,  they  are 
controlled  by  the  Command  Mission  Computer  (CMC)  and  the  Computer  Interface  Unit 
(CIU).  If  RF  energy  is  received  by  any  one  of  the  receivers,  that  information  is  sent  to  the 
Encoder  where  information  is  digitized  and  sent  to  the  CMC  for  further  processing.  Only 
one  receiver’s  data  may  be  processed  at  a  time.  The  system,  therefore,  has  high 
sensitivity,  but  low  probability  of  intercept  with  eight  bands  scanned  by  six  receivers,  and 
only  one  channel  for  processing.  Block  86/89A  aircraft  combine  the  encoder  and  CIU  to 
form  a  CIU/E  which  also  assists  in  matching  emitters  with  known  parameters. 

Additionally,  two  channels  are  provided  to  the  CMC  to  speed  (in  theory,  double) 
processing.  Test  data  does  not  support  increased  processing  speed,  however,  all  aircraft 
will  ultimately  be  deployed  in  the  Block  89A  configuration  and  thus  this  version  will  be  the 
focus  of  this  thesis. 

As  stated  previously,  the  receivers  are  digitally  tuned  and  superhetrodyned. 

Various  methods  are  employed  to  prevent  image  frequencies  from  entering  the  system  for 
further  processing.  However,  these  methods  are  not  fool  proof  and  operators  must  use 
care  to  employ  hardware  fixes,  especially  near  high  impulse  density  energy  sources.  Each 
receiver  has  4-6  antennas  arranged  in  the  tail  of  the  aircraft  to  localize  the  direction  of 
arrival.  For  this  paper  we  will  focus  upon  the  six  antenna  arrangement. 

The  six  antennas  are  arranged  at  60  degree  increments  to  cover  360  degrees  of 
space  around  the  aircraft.  As  signals  are  collected  by  the  receivers,  the  amplitudes  are 
compared  at  the  different  feeds,  and  a  geometric  average  is  sent  to  the  encoder.  This 
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method  of  amplitude  comparison  is  not  very  accurate,  but  in  theory  would  work  as 
follows: 

Assume  at  feed  1  the  amplitude  was  2  units. 

At  feed  2  (clockwise  from  feed  1)  the  amplitude  was  1  unit. 

Therefore,  the  energy  is  twice  as  great  at  feed  1,  so  the  Direction  of  Arrival 
(DOA)  is  closer  to  feed  1,  but  how  much  closer? 

By  taking  a  weighted  average,  the  signal  DOA  must  be  20  degrees  to  tlie 
right  of  feed  1  and  40  degrees  to  the  left  of  feed  2. 

Thus  if  feed  one  is  boresighted  with  the  aircraft  heading,  the  DOA  is  20  degrees  relative 
to  the  aircraft  nose  in  the  direction  of  feed  2  boresight. 

In  theory,  this  system  provides  quick,  accurate  results,  with  inexpensive  hardware. 
Experience,  however,  shows  that  the  six  antenna  system  achieves  a  nominal  +/-10  degree 
accuracy.  Additionally,  a  DOA  does  not  guarantee  sure  location  of  the  detected  emitter. 
Atmospherics  (refraction)  and  multi-path  (reflections  from  ground  obstacles)  contribute 
bogus  DOAs.  Finally,  if  energy  is  detected  at  three  feeds  (which  is  possible  with  close 
proximity  high  powered  signals),  the  detection  is  dropped  and  the  receiver  continues  to  its 
next  programmed  tuning  step. 


D.  EMLOC  FUNCTION 


When  a  receiver  detects  a  threat  (or  friendly  for  that  matter)  emitter  and  a  match  is 
made  with  known  parameters,  it  is  displayed  to  the  operator  in  frequency  and  azimuth 
(after  degrees  true  conversion).  The  operator  may  “hook”  the  displayed  symbol  and  via 
several  software  steps,  select  the  target  emitter  for  Emitter  Location,  or  EMLOC.  The 


OBS  then  runs  each  subsequent  detection  of  this  emitter  match  through  an  algorithm  that 
attempts  to  calculate  the  geographic  location  of  the  emitter. 

The  actual  mathematics  necessary  to  accomplish  this  task  are  not  terribly  difficult 
to  compute,  but  an  explanation  is  necessary.  The  following  chapter  will  describe  the 
geometry  and  basic  assumptions  required  to  compute  a  target  emitter’s  geographic 
location  given  a  change  in  aircraft  position  and  a  change  in  DO  A. 
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n.  THE  PROBLEM 


A.  APPLICATION 

t 

All  DO  As  are  delivered  to  the  OBS  relative  to  the  aircraft.  The  results  need  only 
be  developed  relative  to  the  aircraft  since  navigation  is  now  precise  with  the  advent  of 
Global  Positioning  System  (GPS)  in  Block  89A.  Thus,  with  an  emitter  position  calculated 
relative  to  the  aircraft  position  (which  is  known  with  some  precision)  the  emitter  may  be 
geo-located.  The  triangulation  of  emitters  via  Direction  Finding  (DF)  equipment  may  be 
accomplished  manually,  and  the  “Flat  Earth  Approximation”  geometry  is  developed  quite 
simply  (see  Figure  1). 


Figure  1 :  Flat  Earth  Approximation  Geometry 


where  the  ‘i’  subscript  indicates  the  aircraft  initial  position  at  detection  and  match,  and  the 
‘f  subscript  indicates  the  final  aircraft  position.  The  emitters  coordinates  are  (Xj,  Yt)  are 
unknown.  These  coordinates  can  be  further  transcribed  into  geographic  latitude  and 
longitude,  but  we  will  not  pursue  these  calculations  at  this  time. 

More  importantly,  by  examining  the  flat  earth  approximation  we  will  validate  these 
“simple”  equations,  vice  more  complicated  spherical  solutions  required  by  long  range 
systems.  Since  we  are  only  concerned  with  microwave  frequencies,  we  can  approximate 
the  maximum  distance  for  direct  Electromagnetic  propagation  with: 

d=l23(^  +  ^) 

where  d  =  radar  horizon  in  nautical  miles, 
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h^,h^  =  altitude  in  feet  for  the  aircraft  and  emitter 
assuming  the  4/3  earth  radius  approximation  for  refraction  characteristics.  Thus  if  is 
near  sea  level  and  aircraft  ceiling  is  less  than  45,000  feet,  radar  horizons  should  not  exceed 
261  nautical  miles  (300  statute  miles).  L.Laddie  Coburn  [Ref  1]  proved  that  at  distances 
less  than  360  nautical  miles  the  planar  solution  was  less  than  0.5%  in  error  when 
compared  to  the  exact  spherical  solution,  and  that  for  ranges  less  than  150  nautical  miles, 
the  planar  solution  was  sufficiently  accurate  for  emitter  location  estimation. 

Unfortunately,  the  EMLOC  problem  is  not  quite  this  simple.  With  refractive 
errors  and  multipath,  there  are  many  arriving  DOAs  to  consider.  Also,  the  simple 
geometry  shown  in  Figure  1  is  not  always  so  clear.  Typical  changes  in  bearing  may  only  be 
2-4  degrees  over  several  minutes  flying  time  when  aircraft  aspect  is  nose/tail  and/or  range 
to  target  emitter  is  great.  Therefore,  various  algorithms  are  applied  that  employ  maximum 
likelihood  estimators  to  locate  emitters  and  thus,  most  geographic  locations  are  a 
statistical  approximation  of  the  most  probable  location  of  the  emitter. 

The  algorithm  currently  employed  by  the  EA-6B  is  an  iterative  method  that  was 
developed  by  Daniel  Charles  Opperman  in  June  1977.  He  developed  his  method  and 
refers  to  it  as  the  BOOZOO  algorithm  (not  the  term  most  operators  would  have 
preferred).  During  the  early  1970’s,  several  students  at  the  Naval  Postgraduate  School, 
under  the  tutelage  of  Professor  Harold  Titus  investigated  the  use  of  Kalman  Filters  as  a 
possible  method  for  geo-locating  emitters.  We  will  develop  both  methods  and  evaluate 
their  performance  as  to  accuracy,  processing  speed  and  robustness  under  less  than 
optimum  geometry. 
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ni.  THE  BOOZOO  ALGORITHM 


A.  THE  BOOZOO  ALGORITHM  DEVELOPED 


As  stated  previously,  the  Boozoo  Algorithm  (BZA)  was  developed  by  Opperman 
as  part  of  his  thesis  work  with  the  Tactical  Electronic  Reporting,  Processing,  and 
Evaluation  System  (TERPES)  [Ref  2],  This  system  was  developed  by  the  U.S.  Marine 
Corps  to  take  ES  data  collected  by  a  variety  of  sources  (primarily  the  EA-6B)  and  develop 
the  Electronic  Order  of  Battle  (EOB).  The  system  is  currently  fielded  and  is  organic  to 
the  Marine  Tactical  Electronic  Warfare  Squadrons.  Opperman  developed  his  thesis  based 
upon  equations  provided  by  the  North  American  Rockwell  Corporation  (NAR). 

Opperman  noted  that  the  NAR  equations  were  based  on  a  model  where  the  ES  aircraft 
flew  a  track  long  enough  to  create  at  least  90  degrees  of  change  in  bearing.  This  is  not 
always  possible  so  the  BZA  was  developed  to  account  for  DO  A  error  that  was  significant 
when  compared  to  the  change  in  the  initial  bearing  and  final  bearing.  The  BZA  was 
eventually  incorporated  into  the  TERPES  system  in  the  mid  1980’s  and  into  the  EA-6B 
aircraft  during  the  early  1990’s  with  the  EMLOC  function. 

Opperman  used  the  following  as  a  foundation: 


where  P(a, )  is  the  probability  that  the  ith  aircraft  position  will  register  a  bearing  between 
Oj  +  a,,  and  6^  +  +daj ,  and  is  the  DOA  measurement  error  for  the  zth 
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measurement.  The  value,  cr,  is  the  standard  deviation  of  the  DOA  bearing  error 
associated  with  the  rth  measurement. 

If  we  take  n  independent  measurements,  then  we  have  joint  probability: 


P{a^,...,a„)da^da^...da„  = 


expl 


(-l/2/(£)) 


(2;r)  Jo-,0-,. ..a„ 


-da^d^...da„ 


(3-2) 


where 


(3-3) 


and  (see  figure  2)E  is  the  most  likely  position  of  the  target  emitter,  f  is  the  aircraft 
position  vector,  and  f  is  a  vector  orthogonal  to  f ,  and  lying  in  the  plane  defined  by  the 
vectors,  Q,  (actual  position  vector)  and  E , 
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In  order  to  minimize  the  error,  Opperman  suggested  a  weighting  function  that 
would  be  changed  according  to  where  the  DOA  came  from.  In  our  case  this  would  imply 
some  DOAs  are  more  accurate  than  others.  This  turns  out  to  be  true  in  actuality  and 
therefore  a  weighting  function  is  defined: 


where  k  is  an  arbitrary  constant  of  proportionality. 

Thus  the  probability  would  be  maximized  if  we  minimize 

(3-5) 

As  part  of  his  foundation,  Opperman  does  this,  by  first  making  an  approximation.  If  is 
small  then  sin«.  «  a,  so, 

/(E)  =  k'£lV,Um‘a,  (3-6) 

<=1 

From  Figure  2, 

sin  a,,  w  cos^^  -  ni,  j  /  sin p. 


so  now, 


/[e)  =  k^  (cos'  sin '  p,  j 


(3-7) 


We  may  assume  that 


sin 


fi, 


and  making  another  substitution  with  simplification. 
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(3-8) 


Equation  (3-8)  when  viewed  with  Figure  2  may  be  expressed  in  terms  of  the 
expected  emitter  position,  E . 


(3-9) 


The  function  in  equation  (3-9)  may  be  minimized,  but  some  other  reference  may 
need  to  be  defined.  Assuming  X  and  (j)  are  longitude  and  latitude  of  the  aircraft  position 
respectively,  and  6  remains  the  measured  DOA  as  measured  from  true  north,  we  may 

redefine  the  vectors,  Ejy  and  i  such  that: 

a,  =  cos  A, cos^, 
bf  =  sin  A,  cos^, 
c,  =  sin^^, 

X,  =  sin  A ,  cos  -  sin  ^ ,  cos  A ,  sin 
=  -cos  A,  COS0I  -  sin  sin  A,  sin 
Z,  =  cos^,  sin  9, 


Where.  Finally,  we  use  earth  centered  coordinates  (i,j,k)  to  define: 

E  =  ui  +vj  +  wk 
t,  =  X/  +  Yj  +  Z,k 
a;.  +  bik+c^k 

Now  substituting  into  equation  (3-9) 

f(E)^glh  (3-10) 
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where 


g  =  nk'^fV.^X.u  +  Y.V  +  Z,wf 

J=1 

/z  =  ^  -  c,  v)^  +  {c.u  -  a.\vf  +  {a.v  -  b^uy' 

By  differentiating  with  respect  to  u,  v,  and  w  we  obtain  the  gradient  of  the  function, 
will  have  three  components 

/.=(^. -*,/(£))/*  (3-11) 

/.=(?.-*./(£))/* 

Now,  if  the  gradient  is  allowed  to  approach  zero,  then  will  be  minimized 
and  equations  (3-11)  will  become 

=  Kf[E) 

g.=Kf[E)  (3-12) 

g.  =Kf[E) 

By  differentiating  functions  g  and  h  and  substituting  into  equations  (3-12)  we 


obtain: 


^12 

"l3 

u 

'^1 

^12 

^3' 

u 

^21 

^22 

"23 

V 

=/(£) 

^21 

^22 

^23 

V 

"31 

^32 

«33. 

w 

A. 

^32 

*33. 

w 

where 


(3-13) 
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/=! 

i:=l 


/=1 


1  =  1 


^23  ~  / 


;=1 


and 


*11  =!;(*' +C-) 

i=l 

*12  ~  ~2i^'  *' 

1=1 

/I 

*13  = 


;=1 


*=,  =E(<'?+'r) 

/=1 

*23  ~  “Z*'^' 

1=1 

*,3  =  i(“r+'',’) 


The  solution  for  w,  v,  and  w  is  simplified  if  we  assume  that  the  right  side  of 
equation  (3-13)  has  a  known  value.  Therefore  by  using  a  method  of  iteration  we  obtain: 
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(3-14) 


where  k  -  0, 1,2,3,. ..w,  and  the  zeroth  iteration  of  u,  v,  and  w  are  the  solutions  for  the 
emitters  position. 

Simplifying  the  right  half  of  equation  (3-14), 


«ii 
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(3-15) 


where 


^1  =*ll"*  +*12^  +*13^ 
^2  =  *i2«;t  +*22^'*  +^3^k 

B,  =*,3«t+*23n+*33^t 
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Equation  (3-15)  can  be  solved  by  Cramer’s  Rule,  provided  the  A  matrix  is  singular. 


Then, 
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(3-16) 


This  solution,  however,  leads  to  the  determination  of  the  position  vector  in  both 
direction  and  magnitude.  Since  we  only  need  direction,  we  may  further  simplify  the 
solution. 
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The  remaining  calculations  are  quite  simple  to  find  the  latitude  and  longitude  of  the 
emitter  position  vector,  E  . 
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Longitude  =  arctan(v„,/w„,) 

Latitude  =  arctan 

The  BZA  appears  to  be  a  math  intensive  algorithm  that  would  require  an 
inordinate  amount  of  processing  time  when  compared  to  other  methods  of  estimation.  It 
is  now  appropriate  to  develop  several  alternative  estimators  to  the  same  problem,  the 
Kalman  Filter,  and  Extended  Kalman  Filter. 


16 


IV.  THE  ALTERNATIVES 


A.  THE  KALMAN  FILTER 


The  Kalman  Filter  equations  are  a  mathematical  representation  of  linear  systems 
and  are  adaptive  to  a  variety  of  situations  when  trying  to  simulate,  or  estimate  linear 
systems.  They  have  extension  to  nonlinear  systems,  but  may  not  be  as  accurate  in  certain 
circumstances.  The  application  of  Kalman  Filter  equations  to  the  EMLOC  problem  have 
been  explored  since  the  early  1970’s.  In  June,  1972,  L.  Laddie  Coburn  [Ref  1]  developed 
two  algorithms,  a  Kalman  Filter  and  an  Extended  Kalman  Filter,  since  there  are 
nonlinearities.  Also  of  use  was  the  follow-on  work  to  Coburn  by  Edward  H.  Mills  [Ref 
3]  that  was  useful  in  creating  algorithms  compatible  with  the  AYK-14.  We  will  examine 
both  algorithms. 

The  underlying  mathematical  model  for  the  Kalman  Filter  tracker  consists  of  a 
state  equation  (a  linear  difference  equation)  and  a  linear  measurement  equation.  Noise 
processes  appear  in  both  equations. 

The  state  equation  is: 

A(^  +  l)  =  0(A:  +  p)A(A:)  +  rir(;i:)  (4-1) 

where  +  l|^)is  the  state  transition  matrix,  T  is  a  distribution  matrix  related  to  the 
random  forcing  function,  is  a  random  forcing  function  that  accounts  for  random 
excitations  (noise),  and  k  is  the  time  step. 

The  measurement  equation  is: 
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Z{k)  =  H{k)X{k)  +  V{k)  )t=  1,  2,  3,  ... 


(4-2) 


where  H(k)  is  the  observation  matrix,  and  V{k)  is  the  measurement  noise. 

The  noise  is  assumed  to  have  zero  mean  and  therefore, 
i^(k)vuY]= R(.mkj) 

TI^W(k)W(jYY  =  Q(,k)S(kJ)  (4-3) 

EYi‘WuY]=<^ykj 

where  d(kj)  is  the  Kronecker  delta  function  and  has  value  zero  unless  k  =  j. 

The  Kalman  recursion  equations  are  summarized: 

G{k)  =  P{k\k  -  \)H{ky\H{k)P{k\k  -  \)H{ky  + 

P{k\k)  =  Pik\k  - 1)  -  G{k)H{k)P{k\k  - 1) 

P{k  +  \\k)  =  ^{k  +  \k)P{k\k)^{k  +  \,ky  +0{k)  (4-4) 

Xi^k\k)  =  X{k\k  - 1)  +  G{k)^Z{k)  -  H(,k)X{k\k  - 1)] 

X{k\k  - 1)  =  ^{k,k  -  \)X{k  -  - 1) 

where  X{k\j)  is  the  estimate  of  state  X(/:)  based  on  j  measurements  (Z(l),  Z(2),...Z(j)). 

The  Kalman  filter  gains  are  represented  by  matrix  Gik)  and  P{k\j)  represents  the  error 
covariance  matrix  of  the  estimates. 

Since  we  are  concerned  with  DO  As  and  bearing  rates,  we  will  filter  noisy  observed 
DOAs  to  obtain  an  optimal  estimate.  This  is  best  approximated  by  a  state  transition 
matrix: 


0(A:  +  1,A:) 


1  T{k  + 1) 
0  1 


H(k)  =  [l  0] 


(4-5) 
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Also,  since  DO  As  are  the  only  observed  value,  R(k)  and  W(k)  become  scalars  and  the 
bracketed  terms  in  equations  (4-4)  become  scalars,  also.  This  will  play  a  part  in 
simplifying  and  reducing  the  matrix  equations  to  scalar  form.  The  non-uniformity  of 
arrival  time  for  DO  As  requires  that  the  term  T{k+\)  remain  a  variable  in  each  of  the 
applicable  recursion  equations.  This  requires  that  the  error  covariance  and  gains  be 
computed  continuously  since  they  are  a  function  of  the  sample  interval. 

The  P  and  Q  matrices  have  scalar  components  given  by 


'Qu 

e,r 

.A, 

P22: 

Q2^ 

Qri_ 

and  the  G  matrix  is  a  2  x  1  matrix, 
Rewriting  the  gain  equation  G(k), 


(4-6) 


The  observation  matrix  H(k)  —  [1  0]  forces  the  inverse  term  to  become  a  scalar,  allowing 
us  to  compute  gain  terms  directly: 


(4-7) 


which  results  in  scalar  gain  equations 


G,(k)  = 


G,(k)  = 


P^,(k\k-l)  +  R(k) 
P,,(k\k-l)  +  R(k) 


(4-8) 
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Now  the  first  DOA  is  not  filtered,  but  a  technique  of  filtering  backward  is  available 
to  improve  the  estimate  of  the  initial  DOA.  These  equations  are  similar  to  the  Kalman 
equations,  but  they  update  and  smooth  previous  data  based  on  the  last  DOA  that  has  been 
optimally  estimated.  This  method  was  proposed  by  Rauch  [Ref  4]  in  1963. 

Z(1|A:)  =  1(11^:  -  \)  +  D{\\k)G{k)\Z{k)  -  H{k)X{k\k  -  1)1 

••  (4-9) 

D{\\k)  =  Di\\k  -  \)P{k  -  1|A:  -  l)(D(yl:,A:  -  if  P{k\k  - 1)-' 
where  D(1|A:)  becomes  a  function  of  Tik)  and  must  be  solved  sequentially  before  the 

e{k)^H{k)&{k)+V{k)  (4-10) 

where  D(1|A:)  becomes  a  function  of  T(A:)  and  must  be  solved  sequentially  before  the 
smoothed  estimate  can  be  computed.  These  equations  can  also  be  reduced  to  scalar  form 
and  were  utilized  by  Coburn  in  his  thesis  work. 


Initialization  of  the  filter  requires  knowledge  of  the  system  dynamics,  statistical 
properties  of  the  filter,  and  information  regarding  the  initial  state  of  the  system. 


The  noisy  observation  of  the  state  may  be  made  by 


e{k)  = 


e{k) 


e(k)  = 


e(k) 

e(k) 


where  0(k)  is  a  vector  of  noisy  observations,  and  ©(A)  is  a  state  vector  of  exact  emitter 
bearing  angle  and  bearing  rate.  The  observation  matrix  //(k)  =  [1  0]  is  for  measurement  of 
the  bearing  angle  only.  Bearing  rate  depends  on  speed,  range  and  heading  of  the  aircraft 
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with  respect  to  the  emitter.  Since  the  range  is  initially  unknown,  it  must  be  estimated. 
This  could  be  a  standardized  variable  based  on  the  emitter  to  be  located,  or  could  be  an 
operator  input.  For  this  exercise  it  will  be  denoted  by  the  variable  r.  Therefore, 


6  =  -sin(Relative 
r 

The  remaining  equations  are  simple  substitutions  of  theta  for  X  in  the  basic  Kalman 
equations  presented,  resulting  in  the  optimal  estimation  equation: 


DOA) 


^ 57.29578®/  ^ 
/rad 


V 


3600 


deg/ sec 


(4-11) 


J 


'e[k\k) 

'l 

m' 

6{k  -  \\k  - 1) 

rG,(^)' 

^[k\k) 

0 

1 

/(A:-P-1)J^ 

E(k) 


(4-12) 


where 


E(k)  =  0(k)  -0(k-l\k-\)-  T(k)Sik  -  1|A:  -  1) 


Again,  since  there  is  not  a  uniform  sampling  interval,  the  filter  gains  G  may  vary 
and  not  approach  a  steady  state  value  uniformly. 

Once  the  values  for  theta  are  determined,  an  iterative  approach  to  the  equations  in 
chapter  two  should  lead  to  the  emitter  position. 


B.  THE  EXTENDED  KALMAN  FILTER 


Due  to  the  nonlinearities  of  the  trigonometric  functions,  theoretically  an  Extended 
Kalman  Filter  (EKF)  is  warranted  for  our  situation.  Realistically,  we  will  not  have  a 
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priorii  bearing  rate,  nor  will  our  initial  guesses  at  bearing  rate  be  accurate.  Therefore,  we 
must  face  reality  and  attempt  to  locate  with  only  DOA  data.  If  this  is  the  case,  however,  a 
nonlinear  transformation  must  be  utilized  to  transform  observable  DO  As  into  filtered 
position  estimates. 

This  is  not  as  intimidating  as  it  may  sound,  using  the  Extended  Kalman  Filter 
(EKF)  with  the  help  of  small  perturbation  theory  and  a  Taylor  Series  about  an  initial  point. 
Equation  (4-1)  then  becomes 


Z{k)=N[X{k)-\^V{k)  (4-13) 

where  Z{k)  represents  observable  DO  As,  X{k)  is  the  emitter  position  vector,  and  N 
represents  the  nonlinear  transformation 


6  =  arctan 


(Zj.  -  X)  cos  Lj 


(4-14) 


If  the  position  error  of  the  state  vector  X{k)  is  given  by 


X{k)  =  X{k)-X{_k),  (4-15) 

then  the  true  position  X{k)  may  be  expanded  about  the  most  recent  optimal  estimate, 

X{k)  in  a  Taylor  series  expansion  with  higher  order  terms  thrown  away  as  shown: 


A^[  A^(A:)]  =  A'[l(yt)]  +  M{k)X+. . . 


(4-16) 


where 


M{k)  = 


dX 


for  the  emitter  location  algorithm  given  in  equations  (2-3)  and  (2-4). 
The  Kalman  filter  recursion  equations  become 
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G{k)  =  P{k\k  -  \)M{ky[M{k)P{k\k  -  \)M{kY  + 

P{k\k)  =  P{k\k  - 1)  -  G{k)M{k)P{k\k  - 1)  (4-17) 

P{k  + 1)1  A:)  =  0(A:  +  \,k)P{k\k)^{k  +  1,A:)^  +  Q{k) 

Since  the  emitter  is  almost  always  stationary  (or  near  stationary),  the  matrix  becomes 
the  identity  matrix,  greatly  simplifying  the  equations. 

EKF  filters  can  be  volatile,  and  without  proper  initialization  will  diverge. 

Obviously,  if  the  equations  can  not  be  initialized  a  majority  of  the  time  with  good  results, 
the  EKF  will  be  of  no  use  to  the  operator. 

Now  that  the  various  algorithms,  both  current  and  proposed,  have  been  described 
mathematically,  it  is  time  to  subject  each  to  the  same  scenarios  and  determine  which,  if 
any,  is  clearly  more  reliable,  or  efficient. 
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V.  THE  SIMULATION 


A.  TESTED  PARAMETERS 

The  ultimate  goal  of  this  project  is  to  program  the  three  algorithms  presented  to 
determine  which  method: 

1. )  is  more  computationally  efficient  (fewest  floating  point  operations,  or  flops) 

2. )  requires  least  amount  of  memory  (smallest  matrices  and  arrays  held  in  RAM) 

3. )  converges  to  an  acceptable  range  of  the  target  emitter  in  the  fewest  number  of 
iterations. 

B.  SCENARIO 

A  realistic  simulation  was  developed  to  test  the  algorithms  using  Monte  Carlo 
techniques.  All  computer  code  was  written  in  MATLAB. 

The  simulated  ES  aircraft  flies  North  along  the  Prime  Meridian  (0  deg  Longitude) 
from  1.5  deg  S.  Latitude  to  1.5  deg  N.  Latitude,  a  distance  of  180  Nautical  Miles  (Nmi). 
The  aircraft  will  flies  at  constant  altitude  and  heading  at  a  true  airspeed  of  360  knots. 
Measurements  are  taken  at  12  second  intervals  in  an  attempt  to  geo-locate  the  target 
emitter.  The  true  emitter  location  is  on  the  equator  (0  deg  latitude)  at  1.5  deg  East 
Longitude  Direction  of  arrival  accuracy  was  initially  assumed  to  be  +/-  10  degrees.  A 
representation  of  each  algorithm  is  included  in  the  Appendices. 
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Also,  in  an  effort  to  determine  if  geometry  and/or  sampling  interval  has  an  effect 
on  any  algorithm’s  performance,  the  start  and  stop  points  were  altered. 

The  algorithms  were  tested  as  described  above,  and  a  chapter  is  devoted  to  the 
results  of  each. 
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VI.  BOOZOO  RESULTS 


A.  BOOZOO  ALGORITHM 

The  conversion  of  Opperman’s  code  was  completed  with  the  assistance  of  Mr. 
Lester  Hathaway,  NAWC  Weapons  Division,  NAS  Pt.  Mugu,  CA.  Additionally,  the 
actual  weighting  scheme  for  one  of  the  frequency  bands  was  employed.  Figure  3  is  a 
graphic  representation  of  the  weights  applied  to  each  arriving  signal.  Although  every 
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Figure  3:  EA-6B  Mid-Frequency  Weighting  Scheme 


effort  was  made  to  economize  the  volume  of  computer  code  written  in  support  of  the 
algorithm,  the  AN/A YK- 14  does  not  support  matrix  operations  and  must  therefore  rely 
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upon  scalar  equations.  As  performed  for  this  simulation,  the  BZA  requires  memory 
positions  for  75  variables,  including  37  weights  from  the  weight  look-up  table.  An 
example  of  the  code  utilized  is  included  as  Appendix  A. 


B.  BASIC  SCENARIO 


Figure  4,  depicts  a  typical  ES  track  along  the  Prime  Meridian  and  the  cluster  of 
estimated  positions  for  the  target  emitter  for  one  run  of  the  simulation.  Accuracy  was 
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Figure  4:  Single  Run  of  BOOZOO  Baseline  Simulation 
impressive  as  shown  in  Figure  5  for  both  one  simulation  run,  and  the  average  of  1000 
simulations.  The  number  of  iterative  loops  per  DO  A  cut  was  constant  or  nearly  constant 
at  two  iterative  loops  with  each  loop  comprising  42  floating  point  operations.  Set  up  for 
the  iterative  process  took  an  additional  68  floating  point  operations  for  a  total  of  152 


- 1 - 1 - — 

«  1 

La 

7 

+ 

+ 

+ 

- 1 _ 1 _ 

"f - - 

- — — 1 _ 1 _ 

28 


floating  point  operations  per  typical  DOA  cut.  Finally,  a  simple  time  analysis  was 
accomplished  by  evaluating  the  time  to  complete  one  thousand  iterations  of  the  complete 
scenario.  The  one  thousand  simulations  took  2,256.2  seconds  to  complete,  for  an  average 
of  2.5262  seconds  for  each  loop  through  the  scenario.  Thus,  with  150  measurements 
taken,  the  BZA  process  averages  0.0166  seconds  per  measurement  to  arrive  at  a  position 
estimate.  Figure  5  depicts  the  average  error  associated  with  the  BZA  under  the  conditions 
set  in  Chapter  V.  Figure  6  (shown  on  the  next  page)  is  a  blown  up  depiction  of  average 
error  focusing  upon  the  area  of  convergence  for  the  targeted  emitter. 


Figure  5:  Averaged  Error  1000  Simulations 


The  downward  slope  of  the  average  error  is  very  steep  starting  at  about  50  measurements, 
and  flattens  out  after  75  measurements,  which  corresponds  with  the  aircraft  passing  the 
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abeam  position  to  the  target  emitter.  Measurements  were  terminated  at  150,  which 
corresponds  to  a  theoretical  90  degree  change  in  target  position.  Final  average  range  error 
is  less  than  four  nautical  miles  (Nmi.)  at  the  termination  of  the  test. 


Figure  6:  Averaged  Error  (<20  Nmi)  1000  Simulations 


C.  EFFECTS  OF  SAMPLING  INTERVAL 

A  more  realistic  measurement  interval  in  a  “tactical”  situation  would  be  on  the 
order  of  30-60  seconds  between  measurements.  Therefore,  the  BZA  was  tested  as  in  the 
previous  scenario  and  average  errors  are  depicted  in  Figure  7  and  Figure  8  (page  31)  for 
30  second  interval  and  Figure  9  and  Figure  10  (page  32)  for  the  60  second  interval.  In 
both  cases,  the  steep  slope  is  evident  at  about  one  third  of  total  measurements  with 


30 
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Figure  7;  Averaged  Error  30  Second  Intervals  1000  Loops 


Figure  8:  Averaged  Error(<20  Nmi)  30  Second  Intervals  1000  Loops 
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Error 


Figure  9:  Averaged  Error  60  Second  Intervals  1000  Loops 


Figure  10:  Averaged  Error  (<20  Nmi)  60  Second  Intervals  1000  Loops 


a  leveling  off  following  the  passage  of  the  abeam  position.  Accuracy  at  the  completion  of 
the  scenario  was  under  6  Nmi.  for  the  30  second  interval  and  about  7  Nmi.  for  the  60 
second  interval.  Obviously,  more  data  results  in  more  accuracy,  but  in  the  three  cases 
(Basic,  30  sec.,  60  sec.)  presented,  best  accuracy  is  achieved  just  after  45  degrees  of 
relative  change  in  DO  A. 

D.  EFFECTS  OF  START/STOP  POSITION 

In  an  effort  to  determine  if  the  algorithm  was  affected  by  changes  in  position  for 
start  and  stop  of  the  measurements,  the  basic  scenario  was  modified  so  that  sampling 
interval  remained  12  seconds,  but  the  initial  and  end  points  were  modified.  Figures  1 1 
through  13  (pages  34-35)  show  the  effects  of  starting  at  2  degrees  S.  Lat.  And  continuing 
North  to  the  Equator  for  a  total  of  120  measurements.  Figures  14  through  16  (pages  35- 
36)  depict  the  effects  of  starting  at  the  Equator  and  continuing  North  to  2  degrees  N.  Lat., 
also  for  120  measurements.  The  algorithm  gave  similar  final  results  in  both  cases,  for  an 
average  error  at  the  conclusion  of  the  scenario  of  about  6-7  Nmi.  There  is,  however,  a 
striking  difference  in  the  slope  of  the  range  error  curves,  the  range  error  for  the  Equator 
start  is  much  steeper  and  approaches  near  steady  state  by  the  thirtieth  measurement.  In 
contrast,  the  Equator  finish  scenario  (Figure  12)  requires  over  60  measurements  to 
achieve  the  same  average  accuracy.  Clearly  the  preferred  method  to  locate  the  target  if 
the  90  degree  option  is  not  available  is  to  start  with  the  target  at  least  abeam  the  aircraft 
and  to  continue  on  a  course  perpendicular  to  the  radius  of  an  arc  described  by  the 
theoretical  range  to  the  target.  This  is  due  in  part  to  the  more  rapid  change  in  apparent 
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Figure  11:  2  Deg  S.  -  0  Deg  Single  Run  Estimates 


—  1000X  Simulated  AVG 
*  1X  Simulation 


DOA  Cuts 

Figure  12:  2  Deg  S  -  0  Deg  Averaged  Error  1000  Loops 


Error  (NMi) 


bearing  from  the  aircraft  to  the  target  when  the  target  is  at  the  abeam  position.  Figures  1 1 
and  14  show  the  “memory”  of  the  algorithm  where  it  becomes  apparent  that  the  algorithm 
relies  heavily  upon  the  first  DOA  as  all  estimates  seem  to  come  at  or  near  the  first  relative 
bearing  line. 
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VIL  KALMAN  FILTER  RESULTS 


A.  KALMAN  FILTER 

The  Kalman  Filters  tested  for  this  thesis  are  representations  of  those  developed 
and  already  cited  by  Mills  and  Cobum.  Cobum’s  thesis  relied  on  an  inverse  matrix 
manipulation  that  is  not  supported  by  the  AYK-14.  His  example  is  presented  as  a 
comparison  to  the  actual  tested  example,  as  presented  by  Mills.  Mills  developed  a  scalar 
approximation  to  the  inverse  of  the  P  matrix.  The  same  baseline  scenario  was  mn  with 
both  algorithms  an  yielded  the  results  presented  in  Figure  17. 


0  50  100  150 . 


Figure  17:  Comparison  Between  Mills  and  Cobum  Kalman  Filter 


Thus,  although  Cobum’s  case  is  a  clear  winner,  we  cannot  use  his  example  since  it  cannot 
be  supported  by  the  current  computer  utilized  by  the  typical  ES  aircraft. 


B.  BASIC  SCENARIO 

Figure  18,  depicts  a  typical  ES  track  along  the  Prime  Meridian  and  the  cluster  of 
estimated  positions  for  the  target  emitter  for  one  run  of  the  simulation.  Accuracy  was 


Figure  18:  Single  Run  Kalman  Filter  Baseline  Simulation. 


not  as  impressive  as  the  BZA  and  is  shown  in  Figure  19  for  both  one  simulation  run,  and 
the  average  of  1000  simulations.  The  algorithm  utilizes  1 15  floating  point  operations  per 
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measurement  which  is  almost  a  25  percent  reduction  in  necessary  calculations  when 
compared  to  BZA.  Finally,  a  simple  time  analysis  was  accomplished  by  evaluating  the 
time  to  complete  one  thousand  iterations  of  the  complete  scenario.  The  one  thousand 
simulations  took  1657.2  seconds  to  complete  for  an  average  of  1.6572  seconds  for  each 
loop  through  the  scenario.  Thus,  with  150  measurements  taken,  the  Kalman  Filter  (KF) 
process  averages  0.01 1  seconds  per  measurement  to  arrive  at  a  position  estimate,  which  is 
almost  34  percent  faster  than  the  BZA.  Figure  19  depicts  the  average  error  associated 
with  the  BZA  under  the  conditions  set  in  Chapter  V.  Figure  20  (shown  on  the  next  page) 
is  a  blown  up  depiction  of  average  error  focusing  upon  the  area  of  convergence  for  the 
targeted  emitter. 


Figure  19:  Averaged  Error  1000  Simulations 
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The  downward  slope  of  the  average  error  is  very  steep,  starting  at  about  50  measurements 
and  flattens  out  after  75  measurements,  which  corresponds  with  the  aircraft  passing  the 
abeam  position  to  the  target  emitter.  Measurements  were  terminated  at  150,  which 
corresponds  to  a  theoretical  90  degree  change  in  target  position.  Final  average  range  error 
is  around  25  Nmi.  at  the  termination  of  the  test. 


Figure  20:  Averaged  Error  (<60  Nmi)  1000  Simulations 


C.  EFFECTS  OF  SAMPLING  INTERVAL 


As  in  Chapter  VI,  a  more  realistic  measurement  interval  in  a  “tactical”  situation 
would  be  on  the  order  of  30-60  seconds  between  measurements.  Therefore,  the  Kalman 
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Filter  was  tested  as  in  the  previous  scenario  and  average  errors  are  depicted  in  Figure  21 
and  Figure  22  (page  44)  for  30  second  interval  and  Figure  23  (page  44)  and  Figure  24 
(page  45)  for  the  60  second  interval.  In  both  cases,  the  steep  slope  is  evident  at  about  one 
third  of  total  measurements  with  a  leveling  off  following  the  passage  of  the  abeam 
position.  There  is  a  spike  in  both  scenarios  at  around  the  abeam  position,  but  in  both 
cases  the  algorithm  rapidly  returns  toward  convergence  on  a  minimum  error.  Accuracy  at 
the  completion  of  the  scenario  was  about  26  Nmi.  for  the  30  second  interval  and  about  28 
Nmi.  for  the  60  second  interval.  In  the  case  of  the  Kalman  Filter  it  appears,  more  data 
results  in  more  accuracy,  but  not  a  remarkable  improvement  when  compared  to  the 
baseline.  Again,  in  all  three  cases  (Basic,  30  sec.,  60  sec.)  presented,  best  accuracy  is 
achieved  just  after  45  degrees  of  relative  change  in  DO  A. 


Figure  21:  Averaged  Error  30  Second  Intervals  1000  Loops 


43 


Error  (NMi) 


Figure  23:  Averaged  Error  60  Second  Intervals  1000  Loops 
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Figure  24:  Averaged  Error  (<80  Nmi)  60  Second  Intervals  1000  Loops 


D.  EFFECTS  OF  START/STOP  POSITION 


In  an  effort  to  determine  if  the  Kalman  Filter  algorithm  is  affected  by  changes  in 
position  for  start  and  stop  of  the  measurements,  the  basic  scenario  was  modified  so  that 
sampling  interval  remained  12  seconds,  but  the  initial  and  end  points  were  modified. 
Figures  25  through  27  (pages  45-46)  show  the  effects  of  starting  at  2  degrees  S.  Lat.  And 
continuing  North  to  the  Equator  for  a  total  of  100  measurements.  Figures  28  through  30 
(pages  47-48)  depict  the  effects  of  starting  at  the  Equator  and  continuing  North  to  2 
degrees  N.  Lat.,  also  for  100  measurements.  Both  algorithms  gave  similar  final  results  for 
an  average  error  at  the  conclusion  of  the  scenario  of  about  35  Nmi.  There  is,  however,  a 
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striking  difference  in  the  slope  of  the  range  error  curves,  the  range  error  for  the  Equator 
start  is  much  steeper  and  approaches  near  steady  state  by  the  fiftieth  measurement.  In 
contrast,  the  Equator  finish  scenario  (Figure  26)  requires  over  80  measurements  to 
achieve  the  same  average  accuracy.  Clearly  the  preferred  method  to  locate  the  target,  if 
the  90  degree  option  is  not  available,  is  to  start  with  the  target  abeam  the  aircraft  and  to 
continue  on  a  course  perpendicular  to  the  radius  of  an  arc  described  by  the  theoretical 
range  to  the  target.  This  is  in  agreement  with  results  obtained  utilizing  BZA. 


Figure  25:  2  Deg  S.  -  0  Deg  Single  Run  Estimates 
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Error  (NMi) 


Figure  27:  2  Deg  S  -  0  Deg  Averaged  Error  (<100  NMi)1000  Loops 
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Figure  28;  0  Deg  -  2  Deg  N.  Single  Run  Estimates 


Figure  30:  0  Deg  -  2  Deg  N.  Averaged  Error  (<60  Nmi)  1000  Lx)ops 


Figures  25  and  28  show  the  “memory”  of  the  algorithm  where  it  becomes  apparent  that 
the  algorithm  relies  heavily  upon  the  first  DOA  as  all  estimates  seem  to  come  at  or  near 
the  first  relative  bearing  line.  This  is  expected  with  the  Kalman  Filter  as  the  first  DOA  is 
continuously  smoothed  as  the  algorithm  proceeds,  and  is  consistent  with  BZA. 
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VIIL  EXTENDED  KALMAN  FILTER  RESULTS 

A.  THE  EXTENDED  KALMAN  FILTER 

The  extended  Kalman  Filter  was  programmed  to  mirror  the  work  of  Cobum  and 
Mills,  who  both  relied  upon  a  more  accurate  measuring  device  then  the  one  simulated  for 
this  scenario.  Analysis  of  the  algorithm  was  incomplete  due  to  the  volatility  of  the 
results. 


B.  BASELINE  SCENARIO 

The  results  of  the  baseline  scenario  were  so  poor  that  an  attempt  was  made  to  see 
if  the  problem  was  numerically  unstable  due  to  the  use  of  latitude  and  longitude.  The  use 
of  the  X-Y  plane  could  facilitate  solving  the  equations  and  a  conversion  based  on  aircraft 
present  position  and  heading  could  be  used  to  arrive  at  a  solution.  Unfortunately,  for  DF 
accuracy  greater  than  1  degree,  the  filter  always  diverged.  An  attempt  was  made  to 
provide  only  exact  initialization  data  to  start  the  algorithm,  with  like  results.  Thus,  there 
is  no  apparent  utility  in  this  method  at  this  time. 


51 


52 


IX.  CONCLUSIONS 


Although  the  stated  objective  of  this  thesis  was  to  determine  if  the  Kalman  Filter 
algorithms  would  be  faster  and  more  efficient  than  the  current  implementation,  it  is  noted 
that  accuracy  is  the  ultimate  goal  and  should  not  be  sacrificed  in  order  to  increase  speed, 
and  efficiency. 

A.  EFFICIENCY 

Clearly,  the  evidence  presented  in  chapters  VH  and  Vin  shows  that  the  Kalman 
Filters  are  faster  and  require  fewer  floating  point  operations.  The  AN/A YK- 14  is  a  time 
share  computer  that  has  to  be  optimally  programmed  to  expedite  all  tasks;  however, 
accuracy  remains  paramount. 

B.  MEMORY  REQUIREMENTS 

Again,  the  evidence  shows  the  Kalman  Filter  requires  almost  25  percent  fewer 
memory  registers  in  order  to  complete  the  necessary  equations  to  effectively  locate  target 
emitters.  On  the  other  hand,  the  number  of  required  memory  registers  involved  are 
insubstantial,  and  the  case  can  be  made  that  this  is  not  an  important  criteria  for 
determining  which  algorithm  is  best  utilized. 
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C.  ACCURACY 


Without  doubt,  the  deciding  factor  for  aircraft  operators  given  the  limited 
disparity  in  time  and  memory  requirements  is  effective  accuracy.  The  BZA  is  clearly  far 
more  accurate  than  any  of  the  Kalman  Filters.  Thus,  despite  the  increased  processing 
demands  placed  upon  the  AN/A YK- 14,  the  BZA  should  be  retained.  The  BZA  was  more 
accurate  in  virtually  every  scenario  encountered.  The  effects  of  geometry  and  sample 
interval  were  noticeble,  but  were  still  better  than  both  Kalman  Filters  (Cobum  and  Mills) 
in  the  baseline  scenario. 

D.  FUTURE  DEVELOPMENTS 

As  the  algorithms  present  themselves  now,  the  BZA  is  the  clear  winner.  With  the 
advent  of  portable  terminals  that  interface  with  aircraft  systems,  the  lack  of  speed  and 
computing  power  now  experienced  will  be  overcome  by  the  processing  power  of 
microchips.  It  appears  that  for  this  situation  at  this  time,  the  Kalman  Filter  estimators  are 
not  the  answer.  That  is  not  to  say  a  more  complicated  Kalman  Filter  used  in  conjunction 
with  an  Extended  Kalman  Filter  could  not  be  developed,  but  the  BZA  is  more  accurate 
and  not  remarkably  slower  than  the  currently  programmed  examples  evaluated. 
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E.  AREAS  FOR  FUTURE  RESEARCH 


As  part  of  the  comparison  of  the  various  algorithms,  the  scenarios  and  variables 
have  been  kept  very  simple,  but  true  inflight  experience  has  yielded  results  that  are  almost 
never  as  accurate  as  those  presented  in  this  paper.  I  would  propose  that  the  effect  of 
utilizing  the  BZA  with  an  8-bit  computing  system  inflight  vice  the  32-bit  system  utilized 
for  this  thesis  be  investigated.  Also,  I  have  kept  the  sample  interval  regular  and  not 
extreme  in  length.  If  real  world  constraints  and  limitations  are  considered,  the  results 
may  see  some  impact.  Additionally,  I  have  not  considered  the  various  time  averaging 
program  steps  the  ALQ-99  applies  to  incoming  DO  As  when  the  time  between  reception  is 
long. 

The  EA-6B  is  the  planned  C2W  aircraft  for  the  next  twenty  years,  and  as  such, 
every  effort  must  be  made  to  optimize  the  ES  capability  of  the  crew  as  the  threat  advances 
in  technology  and  processing  power. 
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APPENDIX  A.  BOOZOOCODE 


This  appendix  is  a  representative  case  of  Opperman’s  code  FORTRAN  written  in 
MATLAB. 


%%%%%  BOOZOO.M  %%%%% 


%  Maj  S.P.  Jones,  USMC 

%  This  script  is  a  MATLAB  adaptation  of  the  FORTRAN  subroutine 
%  program  developed  by  Daniel  Opperman  in  June,  1997. 

%  flong  is  aircraft  position  longitude 
%  flat  is  aircraft  position  latitude 

%  azz  is  the  raw  doa  azimuth  bearing  measurement 
%  wt  is  the  weighting  factor  for  the  measurement 
%  i  is  the  index  denoting  the  number  of  the  present  measurement 
%  ncuts  is  the  total  number  of  measurements 
%  elong  is  the  estimated  emitter  longitude  in  radians 
%  elat  is  the  estimated  emitter  latitude  in  radians 


clf 

%%%  Initiate  the  program 

randn(’seed',2)  %  set  seed  generator 

pirad  =  pi/180;  %  conversion  factor  from  deg— >rad 

flong  =  0;  %  0  degrees  long  (rad) 

flat  =  -1.5*pirad;  %  1.5  deg  S.  (rad) 

azzout  =  []; 

flatout=[]; 

flongout  =  []; 

elatmat  =  []; 

elongmat  =  []; 

itout  =  []; 

tgtlatmat  =  []; 

tgtlongmat  =  []; 

range__error  =  []; 

tgtlong  =  1.5*pirad;  %  1.5  deg  E.  (rad) 
tgtlat  =  0;  %  at  the  equator 

wt  =  1; 

wttable  =  [.8  .96  .87  .74  .51  .5  .47  .46  .48  .57  .57  ... 

.62  .42  .64  .61  .64  .62  .67  .67  .65  .6  .62  ... 

.63  .7  .39  .58  .51  .45  .47  .47  .57  .78  .88  ... 

.88  .99  .78  .47]; 
index  =  [0:149]; 

for  i=  1:150 

dferror=  10*pirad*(randn); 
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%azzl  =  atan2((tgtlat-flat),(tgtlong-flong)) 

azz  =  (pi/2)-atan((tgtlat-flat)/(tgtIong-flong))  +  dferror; 

azzout  =  [azzout,azz]; 

%%%  DETERMINE  WEIGHT 

if  azz  <  180 

wtindex  =  round((azz  +  2.5)15); 

else 

wtindex  =  round((360  -  azz  +  2.5)/5); 
end  %%%ifazz<180 

wt  =  wttable(wtindex); 
ifi==l 

test  =  2; 

% 

%%  INITIALIZE  A  AND  B  MATRICES  IF  THIS  IS  FIRST  MEASUREMENT%% 

all  =0; 
al2  =  0; 
al3  =  0; 
a22  =  0; 
a23=0; 
a33  =  0; 

bll  =0; 
bl2  =  0; 
bl3  =  0; 
b22  =  0; 
b23  =  0; 
b33  =  0; 

end  %endifi=l 

%%  CALCULATE  DIRECTIONAL  NUMBERS  XS,YS,ZS  for  Ith  CUTTING  PLANE  %% 

slat  =  sin(flat); 
clat  =  cos(flat); 
slong  =  sin(flong); 
dong  =  cos(flong); 
sazz  =  sin(azz); 
cazz  =  cos(azz); 

XS  =  (slong*cazz  -  slat*dong*sazz)  *  wt; 

YS  =  (-clong*cazz  -  slat*slong*sazz)  *  wt; 

ZS  =  dat*sazz*wt; 

%%%  CALCULATE  THE  A  MATRIX 
all  =all  +  XS^2; 
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al2  =  al2  +  XS*YS; 
al3  =  al3  +  XS*ZS: 
a22  =  a22  +  YS^2; 
a23  =  a23  +  YS*ZS; 
a33  =  a33+ZS''2; 

%  Determine  the  vector  components  X,Y,Z  of  A/C  position  vector 

X  =  clong*clat; 

Y  =  slong  *clat; 

Z  =  slat; 

%  Calculate  the  B  matrix 

bll=bll+Y^2  +  Z^2; 
bl2  =  bl2-X*Y; 
bl3  =  bl3-X*Z; 
b22  =  b22  +  X''2  +  Z''2; 
b23  =  b23  -  Y*Z; 
b33  =  b33+X^2  +  Y^2; 

%%%%  Calculate  the  Co-Factors  of  the  A  Matrix 

cll  =  a22*a33-a23''2; 

Cl2  =  al3*a23-al2*a33; 

Cl3  =  al2*a23-al3*a22; 
c22  =  all*a33-al3''2; 
c23  =  al2*al3-all*a23; 
c33  =  all*a22-al2''2; 

u  =  cll; 
v  =  cl2; 
w  =  cl3; 


elong  =  999; 
elat  =  999; 

if  i  <  2  %%%  Need  at  least  two  cuts 
elat  =  0; 
elong  =  0; 

end  %Ifi<2 

iterations  =  0; 
ifi  =  150 

flops(0) 


end 

while  test  >  1.00000000  %%%Iterative  Loop 
iterations  =  iterations  +1; 


elonl  =  elong; 
elatl  =  elat; 
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A1  =bll*u  +  bl2*v  +  bl3*w; 

A2  =  bl2*u  +  b22*v  +  b23*w; 

A3  =  bl3*u  +  b23*v  +  b33*w; 

%%%  DETERMINE  COMPONENTS  OF  POSITION  VECTOR 

u  =  Al*cl  1  +  A2*cl2  +  A3*cl3; 

V  =  Al*cl2  +  A2*c22  +  A3*c23; 
w  =  Al*cl3  +  A2*c23  +  A3*c33; 

%%%  DETERMINE  LAT  AND  LONG  OF  TGT  EMITTER 

elong  =  atan(v/u); 

elat  =  atan(w/sqrt(u''2  +  v''2)); 

test  =  abs(elonl-elong)  +  abs(elatl-elat); 


end  %%%(while  test  >le-8) 
if  i  ==150 


flops 

end 

itout  =  [itout, iterations]; 
test  =  2; 

fIatout=  [flatout,flat]; 

flongout  =  [flongout,0]; 

flat  =  flat  +  ((6*12/60*l/60)*pirad); 

tgtlatmat  =  [tgtlatinat,tgtlat/pirad]; 
tgtlongmat  =  [tgtlongmat,tgtlong/pirad]; 

elatout  =  elat/pirad; 
elatmat=  [elatmat,elatout]; 

elongout  =  elong/pirad; 
elongmat=  [elongmat,elongout]; 

error  =  sqrt((elat-tgtlat)'^2  +  (elong  -  tgtlong)''2)/pirad*60; 

range_error  =  [range_error, error]; 

end  %%%  for  i  =  1  :  360 

figure(l) 

plot(flongout/pirad,flatout/pirad,'b*’,elongmat,elatmat,'r+',... 

tgtlongmat, flatout/pirad,'y-’,index,tgtlatmat,'y-’) 
xlabel(’Longitude  (deg)  East'),ylabel(’Latitude  (deg)') 
axis([0,3,-1.5,1.5]) 
print  -deps  fig4ch 
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figure(2) 

plot(range_error) 
titleCdistance  error  (Nmi)’) 
print  -deps  bz3fig2 

figure(3) 

plot(elatmat,'+*) 

title('Predicted  Target  Latitude  (deg)  Target  Located  at  0  deg') 
hold 

plot(tgtlatmat) 

axis([0,i,-5,5]) 

figure(4) 

plot(elongmat,'+') 

title('Predicted  Target  Longitude  (deg)  Target  Located  at  1.5  deg') 
hold 

plot(tgtlongmat) 

axis([0,i,~5,5]) 

figure(5) 

plot(itout) 

title('Iterations  Required  to  Make  Estimate') 
axis([0,i,0,max(itout)+2]) 
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APPENDIX  B.  COBURN  KALMAN  FILTER  CODE 


This  appendix  is  a  representative  case  of  Cobum’s  FORTRAN  code  written  in 
MATLAB. 


%%%%%  CKFAVG.M  %%%%% 


%  Maj  S.P.  Jones,  USMC 

%  This  function  is  a  MATLAB  adaptation  of  the  FORTRAN 
%  program  developed  by  L.  Laddie  Cobum  in  June,  1972. 

%  flong  is  aircraft  position  longitude 
%  flat  is  aircraft  position  latitude 
%  vel  is  the  aircraft  velocity 
%  azz  is  the  raw  doa  azimuth  bearing  measurement 
%  wt  is  the  weighting  factor  for  the  measurement(not  used) 

%  i  is  the  index  denoting  the  number  of  the  present  measurement 
%  ncuts  is  the  total  number  of  measurements 
%  toa  is  time  of  arrival 

%  elong  is  the  estimated  emitter  longitude  in  radians 
%  elat  is  the  estimated  emitter  latitude  in  radians 

%  acc  is  df  accuracy  1  sigma  usually  10  degrees 


clear  all 

%%%  Initiate  the  program 
randn(’seed’,2)  %  set  seed  generator 

ncuts  =  150; 
index  =  [l:ncuts-l]; 

degperrad  =  180/pi;  %  conversion  factor  from  deg-“>rad 

elatmatc  =  zeros(l,ncuts-l); 
elongmatc  =  zeros(l, ncuts- 1); 
rngerrsumc  =  zeros(l, ncuts- 1); 

total  =  10; 

for  loops  =  Ltotal 

flatout  =  []; 
flongout  =  []; 
azzout  =  []; 
elatoutc  =  []; 
elongoutc  =  []; 
rngerroroutc  =  []; 
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flongl  =  O/degperrad;  %  0  degrees  long  (rad) 

flatl  =  -1.5/degperrad;  %  1.5  deg  S.  (rad) 
flong  =  flongl; 
flat  =  flatl; 


tgtlong  =  1 .5/degperrad;  %  1 .5  deg  E.  (rad) 

tgtlat  =  O/degperrad;  %  at  the  equator 

tgtlatout  =  tgtlat/degperrad*ones(l,ncuts-l); 
tgtlongout  =1.5*ones(l,ncuts-l); 


acc  =  10/degperrad; 
vel  =  360; 
mgest=  1.414*90; 
toa  =  0; 


%  df  accuracy  is  about  +/- 10  degrees 
%  A/C  velocity  is  360  knots 

%  initial  range  estimate  always  150 


for  i  =  l:ncuts  %travel  180  NMi  measuring  every  12  sec 
dferror  =  acc*(randn); 

azz  =  (pi/2)-atan((tgtlat-flat)/(tgtlong-flong))  +  dferror; 

azzout  =  [azzout,azz]; 
azzdeg  =  azz*degperrad; 

ifi==l 


pll  =  10000; 
pl2  =  0; 
p21=0; 
p22  = 10000; 

qll  =0; 
ql2  =  0; 
q22  =  0; 
olddll  =  1; 
olddl2  =  0; 
oldd21  =  0; 
oldd22=l; 
wll  =  l; 
wl2  =  0; 
w21  =0; 
w22  =  l; 

R  =  (acc)'^2; 


%  P  matrix 


%  Q  matrix 


%  D 


%  (10  degree  std  dev)''2  =  var 


gl  =pll/(pll+R);  %G 
g2  =  pl2/(pll+R); 


thtd  =  azzdeg;  %  thtd  will  be  the  filtered  df  cut 

thtdl  =  azzdeg;  %  thtdl  will  the  smoothed  first  cut 

thetadot  =  vel/mgest*sin(azz)*degperrad/3600;  %  theta  dot  deg/sec 
thetadotl  =  thetadot;  %  smoothed  first  theta  dot 

toa  =  0; 

timel  =  toa;  %  dummy  variable  for  computing  elpased  time 

t  =  0;  %  elapsed  time  between  cuts 


64 


else 


k2  =  0; 
k3  =  l; 


t  =  toa-timel; 

timel  =  toa;  %reset  timel  for  next  cut 

tt  =  t/1000; 

qll  =tt; 

ql2  =  tt''(1.5); 

q22  =  tt''2: 


%%%  KALMAN  RECURSION  EQUATIONS  %%%% 

newpll  =pll*(l-gl)  +  2*pl2*t-  (pl2*gl  +  pll*g2)*t  +  ... 
(p22-pl2*g2)*t^2  +  qll; 

newpl2  =  pl2*(l-gl)  +  (p22  -  pl2*g2)*t  +  ql2; 

newp22  =  p22  -  pl2*g2  +  q22; 

gl  =  newpll/(newpll  +  R); 
g2  =  newpl2/(newpl  1  +  R); 

tptd  =  thtd  +  thetadot*t; 

e  =  azzdeg  -  tptd; 

thtd  =  tptd  +  gl  *e;  %  filtered  df  cut 

thetadot  =  thetadot  +  g2*e;  %  bearing  rate 


%%%  SMOOTHING  EQUATIONS  FOR  FIRST  CUT  %%%% 

%%  Cobum  relied  on  inverse 

%%  matrix  fumction  that  is  not  supported  by  CMS-2. 

P=  [pll  pl2 

p21  p22]; 

Pinv  =  inv(P); 
pinll  =  Pmv(l,l); 
pinl2  =  Pinv(l,2); 
pin22  =  Pinv(2,2); 

qpll  =qll*pinll  +ql2*pinl2; 
qpl2  =  ql  l*pinl2  +  ql2*pin22: 
qp21  =ql2*pinll  +q22*pinl2; 
qp22  =  q  1 2*pin  12  +  q22*pin22; 

dll  =  olddll*(l-qpll)  -  olddl2*qp21  +  olddll*qp21*t; 
dl2  =  olddl2*(l-qp22)  -  olddll*qpl2  +  olddll*(qp22-l)*t; 
d21  =  oldd21*(l-qpll)  -  oldd22*qp21  +  oldd21*qp21*t; 
d22  =  oldd22*(l-qp22)  -  oldd2Pqpl2  +  oldd21*(qp22-l)*t; 
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olddll  =  dll; 
olddl2=:dl2; 
oldd21=d21; 
oldd22  =  d22; 

thtdl  =  thtdl  +  (dll*gl  +  dl2*g2)*e; 
thetadotl  =  thetadotl  +  (d21*gl  +  d22*g2)*e; 

pll  =  newpll; 
pl2  =  newpl2; 
p22  =  newp22; 

thti  =  thtd/degperrad; 
thtli  =  thtdl/degperrad; 

tat  =  tan(thtli)  -  tan  (thti); 
wav  =  (flatl  +  flat)/2; 

tila  =  ((flong'-flongl)*cos(wav)  +  flatl  *tan(thtli)  -  flat*tan(thti))/tat; 
tla  =  ((flong-flongl)*cos(tila)  +  flatl*tan(thtli)  -  flat*tan(thti))/tat; 

elat  =  tla*degperrad; 

tlo  =  flong  +  (tla  -  flat)*tan(thti)/cos(tila); 
elong  =  tlo*degperrad; 


elatoutc  =  [elatoutc,elat]; 
elongoutc  =  [elongoutc, elong]; 

rngerror  =  sqrt((tgtlat-tla)'^2  +  (tgtlong-tlo)'^2)*degperrad*60; 
rngerroroutc  =  [rngerroroutc, rngerror]; 

end  %(ifi==l 

oldlat  =  flat; 
oldlong  =  flong; 
toa  =  toa  +  12; 

flat  =  flat  +  12/60*6/60/degperrad; 
flatout=  [flatout,flat]; 

end  %  for  i  =  1:91 

elatmatc  =  elatmatc  +  elatoutc; 
elongmatc  =  elongmatc  +  elongoutc; 
rngerrsumc  =  rngerrsumc  +  rngerroroutc; 

end  %  for  loops  =  1: 

elatavgc  =  elatmatc/total; 
elongavgc  =  elongmatc/total; 
rngerravgc  =  rngerrsumc/total; 

figure(l) 

plot(index, elatavgc, ’*b', index, tgtlatout/r') 
title('Estimated  Latitude’) 
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axis([0,ncuts,-5,5]) 

figure(2) 

plot(index,elongavgc,'*b',index,tgtIongout,'r') 
title('Estimated  Longitude') 
axis([0,ncuts,0,3]) 

figure(3) 
pIot(rngerravgc) 
titleCError  (Nmi)’) 
axis([0,ncuts,0, 1 00]) 
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APPENDIX  C.  MILLS’  KALMAN  FILTER  CODE 


This  appendix  is  a  representative  case  of  Mills’  FORTRAN  code  written  in 
MATLAB. 


%%%%%  MKFAVG.M  %%%%% 


%  Maj  S.P.  Jones,  USMC 

%  This  function  is  a  MATLAB  adaptation  of  the  FORTRAN 
%  program  developed  by  E.H.  Mills  in  March,  1973. 

%  flong  is  aircraft  position  longitude 
%  flat  is  aircraft  position  latitude 
%  vel  is  the  aircraft  velocity 
%  azz  is  the  raw  doa  azimuth  bearing  measurement 
%  i  denotes  the  number  of  the  present  measurement 
%  toa  is  time  of  arrival 

%  elong  is  the  estimated  emitter  longitude  in  radians 
%  elat  is  the  estimated  emitter  latitude  in  radians 

%  acc  is  df  accuracy  1  sigma  usually  10  degrees 

%  tht  is  angle  theta  True  Azimuth  -  True  Heading  I  have 
%  modeled  the  system  around  an  aircraft  heading  000  and 
%  eliminated  this  calculation 

%  thtl  is  the  smoothed  first  azimuth 

%  thetadot  is  the  bearing  rate  in  rad/sec 
%  thetadotl  is  the  smoothed  first  bearing  rate  in  rad/sec 


clear  all 

%%%  Initiate  the  program 

randnf  seed',2)  %  set  seed  generator 

ncuts  =  150; 
index  =  [l:ncuts-l]; 

degperrad  =  1 80/pi;  %  conversion  factor  from  deg->rad 

elatmat  =  zeros(l,ncuts-l); 
elongmat  =  zeros(l,ncuts-l); 
rngerrsum  =  zeros(l, ncuts- 1); 

total  =  1 ; 

for  loops  =  1  :total 

loops 

flatout  =  []; 
flongout  =  []; 

azzout  =  []; 


69 


elatout  =  []; 
elongout=  []; 
mgerrorout  =  []; 

flongl  =  O/degperrad;  %  0  degrees  long  (rad) 

flatl  =  -1.5/degperrad;  %  1.5  deg  S.  (rad) 
flongs:  flongl; 
flat  =  flatl; 

tgtlong  =  1 .5/degperrad;  %  1 .5  deg  E.  (rad) 
tgtlat  =  O/degperrad;  %  at  the  equator 

tgtlatout  =  tgtlat/degperrad*ones(l,ncuts-l); 
tgtlongout  =1.5*ones(l,ncuts-l); 

acc  =  10/degperrad;  %  df  accuracy  is  about  +/- 10  degrees 

vel  =  360;  %  A/C  velocity  is  360  knots 

mgest  =  1.414*90;  %  initial  range  estimate  always  150 

toa  =  0; 


for  i  =  l:ncuts  %travel  180  NMi  measuring  every  12  sec 
dferror  =  acc*(randn); 

azz  =  (pi/2)-atan((tgtlat-flat)/(tgtlong-flong))  +  dferror; 

azzout  =  [azzout,azz]; 
azzdeg  =  azz*degperrad; 

ifi==l 


pll  =  10000; 
pl2  =  0; 

p22  = 10000; 

qll  =0; 
ql2  =  0; 
q22  =  0; 
dll  =  l; 
wll  =  l; 
wl2  =  0; 
w21=0; 
w22=l; 

R  =  (acc)''2; 


%  P  matrix 


%  Q  matrix 


%  D 


%  (10  degree  std  dev)''2  =  var 


gl  =pll/(pll+R);  %G 
g2  =  pl2/(pll+R); 


tht  =  azz; 
thtl  =  azz; 


thetadot  =  vel/rngest*sin(azz)/3600;  %  theta  dot  deg/sec 

thetadotl  =  thetadot;  %  smoothed  first  theta  dot 

toa  =  0; 


timel  =toa; 
t  =  0; 


%  dummy  variable  for  computing  elpased  time 
%  elapsed  time  between  cuts 
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else 


t  =  toa-  timel; 

timel  =  toa;  %reset  timel  for  next  cut 

tt  =  t/1000; 

qll  =tt; 

ql2  =  ttM.5); 

q22  =  tt^2; 


%%%  KALMAN  RECURSION  EQUATIONS  %%%% 

newpll  =  pll*(l-gl)  +  2*pl2*t-  (pl2*gl  +  pll*g2)*t  +  ... 
(p22-pl2*g2)*t^2  +  qll; 

newpl2  =  pl2*(l“gl)  +  (p22  -  pl2*g2)*t  +  ql2; 

newp22  =  p22  -  pl2*g2  +  q22; 

newgl  newpl  l/(newpl  1  +  R); 
newg2  =  newpl2/(newpll  +  R); 

tpt  =  tht  +  thetadot*t; 


e  =  azz  -  tpt; 


tht  =  tpt  +  gl *e;  %  filtered  df  cut 

thetadot  =  thetadot  +  g2*e;  %  bearing  rate 


%%%  SMOOTHING  EQUATIONS  FOR  FIRST  CUT  %%%% 
%%  These  equations  are  extracted  from  Edward  H.  Mills  thesis 
%%  that  builds  upon  Coburn's  work.  Cobum  relied  on  inverse 
%%  matrix  fumction  that  is  not  supported  by  CMS-2,  and  so  I 
%%  have  used  this  group  of  equations  as  developed  by  Mills 
%%  to  better  simulate  operation  onboard  an  actual  aircraft. 

smoothl  =  1  -  pll*(l-gl)/R; 
smooth2  =  smooth  l*t; 
smooths  =  -pl2  *  (l-gl)/R; 
smooth4  =  1  +  smooth3*t; 

newwll  =  wll*smoothl  +  wl2*smooth2; 
newwl2  =  wl  l*smooth3  +  wl2*smooth4; 
neww21  =  w21*smoothl  +  w22*smooth2; 
neww22  =  w21*smooth3  +  w22*smooth4; 

wll  =  newwll; 
wl2  =  newwl2; 
w21  =  neww21; 
w22  =  neww22; 

dll  =  dll  -  (wll''2)*(pll  +R)/(R''2); 
thtl  =thtl  +  (wll/R)*e; 
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thetadotl  =  thetadotl  +  (w21/R)*e; 


pll  =  newpll; 
pl2  =  newpl2; 
p22  =  newp22; 

gl  =newgl; 
g2  =  newg2; 


thti  =  tht; 
thtli  =  thtl; 

tat  =  tan(thtli)  -  tan  (thti); 
wav  =  (flatl  +  flat)/2; 

tila  =  ((flong-flongl)*cos(wav)  +  flatl  *tan(thtli)  -  flat*tan(thti))/tat; 
tla  =  ((flong-flongl)*cos(tila)  +  flatl*tan(thtli)  -  flat*tan(thti))/tat; 
tlo  =  flong  +  (tla  -  flat)*tan(thti)/cos(tila); 

elat  =  tla*degperrad; 

elong  =  tlo*degperrad; 

elatout  =  [elatout,elat]; 
elongout  =  [elongout, elong]; 

rngerror  =  (sqrt((tgtlat-tla)'^2  +  (tgtlong-tlo)''2))*degperrad*60; 
mgerrorout  =  [mgerrorout,mgerror]; 

end  %(if  i=l 

oldlat  =  flat; 
oldlong  =  flong; 
toa  =  toa  +12; 

flat  =  flat  +  12/60*6/60/degperrad; 
flatout  =  [flatout,flat]; 

end  %  for  i  =  l:ncuts 

elatmat  =  elatmat  +  elatout; 
elongmat  =  elongmat  +  elongout; 
rngerrsum  =  rngeirsum  +  mgerrorout; 

end  %  for  loops  =  1: 

elatavg  =  elatmat/total; 
elongavg  =  elongmat/total; 
mgerravg  =  rngerrsum/total; 

figure(l) 

plot(index, elatavg, '*b',index,tgtlatout,’r') 

%title('Estimated  Latitude  150  Cuts  10+/-Degree  Accuracy  Mills  Kalman  Filter') 

xlabel('Cuts'),  ylabel('Latitude  (Deg)') 

axis([0,ncuts,-5,5]) 

%print  -deps  inkfl2fl 


figure(2) 

plot(index,elongavg,'*b',index,tgtlongout,'r') 

%title('Estimated  Longitude  150  Cuts  10+/-Degree  Accuracy  Mills  Kalman  Filter') 

xlabel('Cuts'),  ylabel('Longitude  (Deg)') 

axis([0,ncuts,0,3]) 

%print  -deps  nikfl2f2 

figure(3) 

plot(index,mgerravg, ’-'.index, mgerrorout,'*') 
legend('1000X  Simulated  AVG'.’IX  Simulation') 

%title('Error  (NMi)  150  Cuts  10+/-Degree  Accuracy  Mills  Kalman  Filter') 
xlabel('Cuts'),ylabel('Error  (NMi)') 
axis([0,ncuts, 0,200]) 

%print  -deps  mkfl2f3 

figure(4) 

plot(elongout,elatout,'*g’) 

ylabel('Latitude  (Deg)'),  xlabel('Longitude  (Deg)') 

grid 

%title(’Mills  Kalman  Filter  Predictions  150  Cuts  +/- 10  Degree  Accuracy') 
axis([0,3,-.5,1.5]) 

%print  -deps  mkfl2f4 

figure(5) 

plot(index,mgerravg,’-',index,rngerrorout,'*') 
legendClOOOX  Simulated  AVG’.’IX  Simulation’) 

%title('Error  (NMi)  150  Cuts  10+/-Degree  Accuracy  Mills  Kalman  Filter') 

axis([0,ncuts,0,60]) 

xlabel(’Cuts'),ylabel('Error  (NMi)') 

%print  -deps  mkf  1 2f5 
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